TeaTime4schools: Joint analysis - bacteria
07 Beta diversity analysis
Roey Angel
2021-05-21
Beta diversity analysis
This analysis explores the alpha-diversity ditribution patters in the different samples, based on the DADA2-produced sequences.
Setting general parameters:
set.seed(1000)
min_lib_size <- 5000
data_path <- "./DADA2_pseudo/"
Ps_file <- "TeaTime4Schools_filt3_wTree.Rds"
Proj_name <- "TeaTime4Schools"Load phyloseq object
This phyloseq object was created in 06_Phylogentic_analysis.html. The Ps_obj_filt object excludes contaminants and all sequences classified as eukaryota, chloroplast, mitochondria or unknown but still includes taxa with low prevalence
readRDS(file = paste0(data_path, Ps_file)) %>%
subset_samples(., Field != "Unburied") %>% # drop unburied samples
prune_samples(sample_sums(.) > min_lib_size, .) %>% # remove samples < min_lib_size
filter_taxa(., function(x) sum(x) > 0, TRUE) -> # drop taxa with 0 abundance
Ps_obj_filtqplot(rowSums(otu_table(Ps_obj_filt)), geom = "histogram") +
xlab("Library size")
qplot(log10(rowSums(otu_table(Ps_obj_filt)))) +
xlab("Logged library size")
Standardize abundances to the median sequencing depth (and convert to proportion)
adonis(
otu_table(Ps_obj_filt) ~ Lib.size,
data =
as(sample_data(Ps_obj_filt), "data.frame"),
method = "horn",
permutations = 999
)##
## Call:
## adonis(formula = otu_table(Ps_obj_filt) ~ Lib.size, data = as(sample_data(Ps_obj_filt), "data.frame"), permutations = 999, method = "horn")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Lib.size 1 6.080 6.0804 19.569 0.14122 0.001 ***
## Residuals 119 36.975 0.3107 0.85878
## Total 120 43.056 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ps_obj_filt %>%
otu_table(.) %>%
as(., "matrix") %>%
rowSums() %>%
median() ->
total
standf = function(x, t = total) round(t * (x / sum(x)))
Ps_obj_filt_median <- transform_sample_counts(Ps_obj_filt, standf) # Standardize abundances to median sequencing depth
Ps_obj_filt_median_rel <- transform_sample_counts(Ps_obj_filt_median, function(x) x / sum(x)) # convert to relative abundance (just incase it's explicitly needed)
sample_data(Ps_obj_filt_median)$Lib.size <- sample_sums(Ps_obj_filt_median)
qplot(rowSums(otu_table(Ps_obj_filt_median)), geom = "histogram") +
xlab("Library size")
PlotLibDist(Ps_obj_filt_median)
PlotReadHist(as(otu_table(Ps_obj_filt_median), "matrix"))
notAllZero <- (rowSums(t(otu_table(Ps_obj_filt_median))) > 0)
meanSdPlot(as(log2(t(otu_table(Ps_obj_filt_median))[notAllZero, ] + 1), "matrix"))
Variance partitioning models and ordinations
Partitioning the data using discrete distance
adonis(
otu_table(Ps_obj_filt_median) ~ Lib.size,
data =
as(sample_data(Ps_obj_filt_median), "data.frame"),
method = "horn",
permutations = 999
)##
## Call:
## adonis(formula = otu_table(Ps_obj_filt_median) ~ Lib.size, data = as(sample_data(Ps_obj_filt_median), "data.frame"), permutations = 999, method = "horn")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Lib.size 1 0.285 0.28485 0.79255 0.00662 0.57
## Residuals 119 42.770 0.35941 0.99338
## Total 120 43.055 1.00000adonis(
otu_table(Ps_obj_filt_median_rel) ~ Field + Sample.type * Season,
data =
as(sample_data(Ps_obj_filt_median_rel), "data.frame"),
method = "horn",
permutations = 999
)##
## Call:
## adonis(formula = otu_table(Ps_obj_filt_median_rel) ~ Field + Sample.type * Season, data = as(sample_data(Ps_obj_filt_median_rel), "data.frame"), permutations = 999, method = "horn")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Field 2 1.158 0.5789 5.305 0.02689 0.001 ***
## Sample.type 2 16.495 8.2477 75.584 0.38313 0.001 ***
## Season 3 8.674 2.8913 26.497 0.20146 0.001 ***
## Sample.type:Season 6 5.052 0.8419 7.716 0.11733 0.001 ***
## Residuals 107 11.676 0.1091 0.27119
## Total 120 43.055 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# sample_data(Ps_obj_filt_median)$Type.season <- paste(sample_data(Ps_obj_filt)$Sample.type, sample_data(Ps_obj_filt)$Season) # too many comparisons
# Ps_obj_filt_median_s <- prune_taxa(names(sort(taxa_sums(Ps_obj_filt_median), TRUE)[1:100]), Ps_obj_filt_median) # for testing
# Compare sample-type pairs
mod_pairwise1 <- PairwiseAdonis(
otu_table(Ps_obj_filt_median),
sample_data(Ps_obj_filt_median)$Sample.type,
sim.function = "vegdist",
sim.method = "horn",
p.adjust.m = "BH"
)
print(mod_pairwise1)## pairs total.DF F.Model R2 p.value p.adjusted sig
## 1 Green tea vs Rooibos 91 13.90875 0.1338555 0.001 0.001 **
## 2 Green tea vs Soil 73 48.04850 0.4002424 0.001 0.001 **
## 3 Rooibos vs Soil 75 67.00992 0.4752142 0.001 0.001 **(sig_pairs1 <- as.character(mod_pairwise1$pairs[mod_pairwise1$p.adjusted < 0.05]))## [1] "Green tea vs Rooibos" "Green tea vs Soil" "Rooibos vs Soil"# (meaningful_sig_pairs1 <- c("Green tea vs Rooibos", "Green tea vs Soil", "Rooibos vs Soil"))
# compare seasons
mod_pairwise2 <- PairwiseAdonis(
otu_table(Ps_obj_filt_median),
sample_data(Ps_obj_filt_median)$Season,
sim.function = "vegdist",
sim.method = "horn",
p.adjust.m = "BH"
)
print(mod_pairwise2)## pairs total.DF F.Model R2 p.value p.adjusted sig
## 1 Winter vs Spring 58 7.531288 0.11670754 0.001 0.0012 *
## 2 Winter vs Summer 59 21.038047 0.26617620 0.001 0.0012 *
## 3 Winter vs Autumn 61 16.210974 0.21271181 0.001 0.0012 *
## 4 Spring vs Summer 58 7.851482 0.12106866 0.001 0.0012 *
## 5 Spring vs Autumn 60 6.103328 0.09374832 0.001 0.0012 *
## 6 Summer vs Autumn 61 5.653014 0.08610441 0.002 0.0020 *(sig_pairs2 <- as.character(mod_pairwise2$pairs[mod_pairwise2$p.adjusted < 0.05]))## [1] "Winter vs Spring" "Winter vs Summer" "Winter vs Autumn" "Spring vs Summer"
## [5] "Spring vs Autumn" "Summer vs Autumn"# (meaningful_sig_pairs2 <- c("Green tea vs Rooibos", "Green tea vs Soil", "Rooibos vs Soil"))Calculate ordinations
Ps_obj_ord1 <- ordinate(Ps_obj_filt_median, "CAP", "horn", formula = Ps_obj_filt_median ~ Field + Sample.type * Season)
evals <- eigenvals(Ps_obj_ord1) # /sum(eigenvals(Ps_obj_ord)) * 100
Ps_obj_filt_median %>%
plot_ordination(., Ps_obj_ord1, type = "samples", shape = "Field", color = "Sample.type", justDF = TRUE) %>%
mutate_at(., "Season", ~fct_relevel(., "Winter", "Spring", "Summer", "Autumn")) %>%
mutate_at(., "Sample.type", ~fct_relevel(., "Soil", "Green tea", "Rooibos")) %>%
dplyr::rename(., `Sample type` = Sample.type) ->
ord_df
# p_ord.file <- "PCoA_bray"
# svglite(paste0(p_ord.file, ".svg"),
# width = 10, height = 7.2)
p_ord <- ggplot(ord_df,
aes(
x = CAP1,
y = CAP2,
shape = Field,
color = `Sample type`
)) +
geom_point(size = 4, alpha = 2 / 3) +
theme_bw(base_size = 14) +
scale_colour_manual(values = pal_locuszoom("default")(3)[c(2, 3, 1)]) +
stat_ellipse(
aes(x = CAP1, y = CAP2, group = `Sample type`),
alpha = 0.5,
type = "norm",
level = 0.95,
linetype = 2,
inherit.aes = FALSE
) +
labs(x = sprintf("CAP1 [%s%%]", round(evals[1], 1)),
y = sprintf("CAP2 [%s%%]", round(evals[2], 1))) +
coord_fixed(sqrt(evals[2] / evals[1])) +
facet_wrap(~ Season)
# # Now add the environmental variables as arrows
# arrowmat <- scores(Ps_obj_ord1, display = "bp") # bipplot arrows
# # Add labels, make a data.frame
# arrowdf <- data.frame(labels = rownames(arrowmat), arrowmat)
#
# arrowdf %<>%
# mutate(labels = fct_recode(labels,
# "Braunerde" = "FieldBraunerde",
# "Kolluvisol" = "FieldKolluvisol",
# "Rooibos" = "Sample.typeRooibos",
# "Soil" = "Sample.typeSoil"
# ))
# arrowdf %<>% dplyr::slice(., c(1:4))
#
# # Define the arrow aesthetic mapping
# arrow_map = aes(xend = CAP1, yend = CAP2, x = 0, y = 0, shape = NULL, color = NULL,
# label = labels)
# label_map = aes(x = 1.2 * CAP1, y = 1.2 * CAP2, shape = NULL, color = NULL,
# label = labels)
# # Make a new graphic
# arrowhead = arrow(length = unit(0.05, "npc"))
# p_ord <- p_ord +
# geom_segment(
# arrow_map,
# size = 0.5,
# data = arrowdf,
# color = "gray",
# arrow = arrowhead
# ) +
# geom_text(label_map, size = 2, data = arrowdf)
print(p_ord)
# dev.off()
# ggsave(
# paste0(p_ord.file, ".png"),
# p_ord,
# device = "png",
# width = 10,
# height = 6
# )
# gz(paste0(p_ord.file, ".svg"), paste0(p_ord.file, ".svgz"))Partitioning the data using phylometric distance
Unifrac_mat <- UniFrac(Ps_obj_filt_median,
weighted = TRUE,
normalized = TRUE,
parallel = TRUE,
fast = TRUE)
adonis(
Unifrac_mat ~ Lib.size,
data =
as(sample_data(Ps_obj_filt_median), "data.frame"),
permutations = 999
)##
## Call:
## adonis(formula = Unifrac_mat ~ Lib.size, data = as(sample_data(Ps_obj_filt_median), "data.frame"), permutations = 999)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Lib.size 1 0.0589 0.058917 1.2597 0.01047 0.254
## Residuals 119 5.5658 0.046771 0.98953
## Total 120 5.6247 1.00000adonis(
Unifrac_mat ~ Field + Sample.type * Season,
data =
as(sample_data(Ps_obj_filt_median), "data.frame"),
permutations = 999
)##
## Call:
## adonis(formula = Unifrac_mat ~ Field + Sample.type * Season, data = as(sample_data(Ps_obj_filt_median), "data.frame"), permutations = 999)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Field 2 0.0703 0.03515 3.295 0.01250 0.007 **
## Sample.type 2 2.6715 1.33576 125.228 0.47496 0.001 ***
## Season 3 1.1516 0.38388 35.989 0.20475 0.001 ***
## Sample.type:Season 6 0.5899 0.09832 9.218 0.10488 0.001 ***
## Residuals 107 1.1413 0.01067 0.20291
## Total 120 5.6247 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Calculate ordinations
Ps_obj_ord2 <- ordinate(Ps_obj_filt_median, "PCoA", "unifrac", weighted = TRUE, formula = Ps_obj_filt_median ~ Field + Sample.type * Season)
evals <- Ps_obj_ord2$values$Eigenvalues[1:2] # /sum(eigenvals(Ps_obj_ord)) * 100
Ps_obj_filt_median %>%
plot_ordination(., Ps_obj_ord2, type = "samples", shape = "Field", color = "Sample.type", justDF = TRUE) %>%
mutate_at(., "Season", ~fct_relevel(., "Winter", "Spring", "Summer", "Autumn")) %>%
mutate_at(., "Sample.type", ~fct_relevel(., "Soil", "Green tea", "Rooibos")) %>%
dplyr::rename(., `Sample type` = Sample.type) ->
ord_df2
# p_ord2.file <- "PCoA_bray"
# svglite(paste0(p_ord2.file, ".svg"),
# width = 10, height = 7.2)
p_ord2 <- ggplot(ord_df2,
aes(
x = Axis.1,
y = Axis.2,
shape = Field,
color = `Sample type`
)) +
geom_point(size = 4, alpha = 2 / 3) +
theme_bw(base_size = 14) +
scale_colour_manual(values = pal_locuszoom("default")(3)[c(2, 3, 1)]) +
stat_ellipse(
aes(x = Axis.1, y = Axis.2, group = `Sample type`),
alpha = 0.5,
type = "norm",
level = 0.95,
linetype = 2,
inherit.aes = FALSE
) +
labs(x = sprintf("Axis1 [%s%%]", round(evals[1], 1)),
y = sprintf("Axis2 [%s%%]", round(evals[2], 1))) +
coord_fixed(sqrt(evals[2] / evals[1])) +
facet_wrap( ~ Season)
# plot_grid(p_ord,p_ord2)
print(p_ord2)
# dev.off()
# ggsave(
# paste0(p_ord2.file, ".png"),
# p_ord2,
# device = "png",
# width = 10,
# height = 6
# )
# gz(paste0(p_ord2.file, ".svg"), paste0(p_ord2.file, ".svgz"))Variance partitioning models and ordinations - soils excluded
Partitioning the data using discrete distance
Ps_obj_filt_median_tea <- subset_samples(Ps_obj_filt_median, Sample.type != "Soil")
adonis(
otu_table(Ps_obj_filt_median_tea) ~ Lib.size,
data =
as(sample_data(Ps_obj_filt_median_tea), "data.frame"),
method = "horn",
permutations = 999
)##
## Call:
## adonis(formula = otu_table(Ps_obj_filt_median_tea) ~ Lib.size, data = as(sample_data(Ps_obj_filt_median_tea), "data.frame"), permutations = 999, method = "horn")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Lib.size 1 0.3309 0.33093 1.0424 0.01145 0.416
## Residuals 90 28.5712 0.31746 0.98855
## Total 91 28.9021 1.00000adonis(
otu_table(Ps_obj_filt_median_tea) ~ Field + Sample.type * Season,
data =
as(sample_data(Ps_obj_filt_median_tea), "data.frame"),
method = "horn",
permutations = 999
)##
## Call:
## adonis(formula = otu_table(Ps_obj_filt_median_tea) ~ Field + Sample.type * Season, data = as(sample_data(Ps_obj_filt_median_tea), "data.frame"), permutations = 999, method = "horn")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Field 2 1.3803 0.6901 5.4738 0.04776 0.001 ***
## Sample.type 1 3.8840 3.8840 30.8053 0.13439 0.001 ***
## Season 3 10.4768 3.4923 27.6981 0.36249 0.001 ***
## Sample.type:Season 3 2.8222 0.9407 7.4612 0.09765 0.001 ***
## Residuals 82 10.3388 0.1261 0.35772
## Total 91 28.9021 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# sample_data(Ps_obj_filt_median)$Type.season <- paste(sample_data(Ps_obj_filt)$Sample.type, sample_data(Ps_obj_filt)$Season) # too many comparisons
# Ps_obj_filt_median_s <- prune_taxa(names(sort(taxa_sums(Ps_obj_filt_median), TRUE)[1:100]), Ps_obj_filt_median) # for testing
# Compare sample-type pairs
mod_pairwise3 <- PairwiseAdonis(
otu_table(Ps_obj_filt_median_tea),
sample_data(Ps_obj_filt_median_tea)$Sample.type,
sim.function = "vegdist",
sim.method = "horn",
p.adjust.m = "BH"
)
print(mod_pairwise3)## pairs total.DF F.Model R2 p.value p.adjusted sig
## 1 Green tea vs Rooibos 91 13.90875 0.1338555 0.001 0.001 **(sig_pairs3 <- as.character(mod_pairwise3$pairs[mod_pairwise3$p.adjusted < 0.05]))## [1] "Green tea vs Rooibos"# (meaningful_sig_pairs3 <- c("Green tea vs Rooibos", "Green tea vs Soil", "Rooibos vs Soil"))
# compare seasons
mod_pairwise4 <- PairwiseAdonis(
otu_table(Ps_obj_filt_median_tea),
sample_data(Ps_obj_filt_median_tea)$Season,
sim.function = "vegdist",
sim.method = "horn",
p.adjust.m = "BH"
)
print(mod_pairwise4)## pairs total.DF F.Model R2 p.value p.adjusted sig
## 1 Winter vs Spring 46 13.19236 0.2267026 0.001 0.001 **
## 2 Winter vs Summer 44 34.98156 0.4485876 0.001 0.001 **
## 3 Winter vs Autumn 47 28.36909 0.3814635 0.001 0.001 **
## 4 Spring vs Summer 43 11.80495 0.2194027 0.001 0.001 **
## 5 Spring vs Autumn 46 10.05915 0.1826972 0.001 0.001 **
## 6 Summer vs Autumn 44 8.91082 0.1716563 0.001 0.001 **(sig_pairs4 <- as.character(mod_pairwise4$pairs[mod_pairwise4$p.adjusted < 0.05]))## [1] "Winter vs Spring" "Winter vs Summer" "Winter vs Autumn" "Spring vs Summer"
## [5] "Spring vs Autumn" "Summer vs Autumn"# (meaningful_sig_pairs2 <- c("Green tea vs Rooibos", "Green tea vs Soil", "Rooibos vs Soil"))Calculate ordinations
Ps_obj_ord3 <- ordinate(Ps_obj_filt_median_tea, "CAP", "horn", formula = Ps_obj_filt_median_tea ~ Field + Sample.type * Season)
evals <- eigenvals(Ps_obj_ord3) # /sum(eigenvals(Ps_obj_ord)) * 100
Ps_obj_filt_median_tea %>%
plot_ordination(., Ps_obj_ord3, type = "samples", shape = "Field", color = "Sample.type", justDF = TRUE) %>%
mutate_at(., "Season", ~fct_relevel(., "Winter", "Spring", "Summer", "Autumn")) %>%
mutate_at(., "Sample.type", ~fct_relevel(., "Green tea", "Rooibos")) %>%
dplyr::rename(., `Sample type` = Sample.type) ->
ord_df3
# p_ord.file <- "PCoA_bray"
# svglite(paste0(p_ord.file, ".svg"),
# width = 10, height = 7.2)
p_ord3 <- ggplot(ord_df3,
aes(
x = CAP1,
y = CAP2,
shape = Field,
color = `Sample type`
)) +
geom_point(size = 4, alpha = 2 / 3) +
theme_bw(base_size = 14) +
scale_colour_manual(values = pal_locuszoom("default")(3)[c(3, 1)]) +
stat_ellipse(
aes(x = CAP1, y = CAP2, group = `Sample type`),
alpha = 0.5,
type = "norm",
level = 0.95,
linetype = 2,
inherit.aes = FALSE
) +
labs(x = sprintf("Axis1 [%s%%]", round(evals[1], 1)),
y = sprintf("Axis2 [%s%%]", round(evals[2], 1))) +
coord_fixed(sqrt(evals[2] / evals[1])) +
facet_wrap(~ Season)
print(p_ord3)
# dev.off()
# ggsave(
# paste0(p_ord.file, ".png"),
# p_ord,
# device = "png",
# width = 10,
# height = 6
# )
# gz(paste0(p_ord.file, ".svg"), paste0(p_ord.file, ".svgz"))Partitioning the data using phylometric distance
Unifrac_mat <- UniFrac(Ps_obj_filt_median_tea,
weighted = TRUE,
normalized = TRUE,
parallel = TRUE,
fast = TRUE)
adonis(
Unifrac_mat ~ Lib.size,
data =
as(sample_data(Ps_obj_filt_median_tea), "data.frame"),
permutations = 999
)##
## Call:
## adonis(formula = Unifrac_mat ~ Lib.size, data = as(sample_data(Ps_obj_filt_median_tea), "data.frame"), permutations = 999)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Lib.size 1 0.0431 0.043104 2.0314 0.02207 0.076 .
## Residuals 90 1.9097 0.021219 0.97793
## Total 91 1.9528 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(
Unifrac_mat ~ Field + Sample.type * Season,
data =
as(sample_data(Ps_obj_filt_median_tea), "data.frame"),
permutations = 999
)##
## Call:
## adonis(formula = Unifrac_mat ~ Field + Sample.type * Season, data = as(sample_data(Ps_obj_filt_median_tea), "data.frame"), permutations = 999)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Field 2 0.04872 0.024358 3.170 0.02495 0.008 **
## Sample.type 1 0.16444 0.164441 21.403 0.08421 0.001 ***
## Season 3 0.94795 0.315983 41.127 0.48543 0.001 ***
## Sample.type:Season 3 0.16168 0.053893 7.015 0.08279 0.001 ***
## Residuals 82 0.63001 0.007683 0.32262
## Total 91 1.95279 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Calculate ordinations
Ps_obj_ord4 <- ordinate(Ps_obj_filt_median_tea, "PCoA", "unifrac", weighted = TRUE, formula = Ps_obj_filt_median_tea ~ Field + Sample.type * Season)
evals <- Ps_obj_ord4$values$Eigenvalues[1:2] # /sum(eigenvals(Ps_obj_ord)) * 100
Ps_obj_filt_median_tea %>%
plot_ordination(., Ps_obj_ord4, type = "samples", shape = "Field", color = "Sample.type", justDF = TRUE) %>%
mutate_at(., "Season", ~fct_relevel(., "Winter", "Spring", "Summer", "Autumn")) %>%
mutate_at(., "Sample.type", ~fct_relevel(., "Green tea", "Rooibos")) %>%
dplyr::rename(., `Sample type` = Sample.type) ->
ord_df4
# p_ord.file <- "PCoA_bray"
# svglite(paste0(p_ord.file, ".svg"),
# width = 10, height = 7.2)
p_ord4 <- ggplot(ord_df4,
aes(
x = Axis.1,
y = Axis.2,
shape = Field,
color = `Sample type`
)) +
geom_point(size = 4, alpha = 2 / 3) +
theme_bw(base_size = 14) +
scale_colour_manual(values = pal_locuszoom("default")(3)[c(3, 1)]) +
stat_ellipse(
aes(x = Axis.1, y = Axis.2, group = `Sample type`),
alpha = 0.5,
type = "norm",
level = 0.95,
linetype = 2,
inherit.aes = FALSE
) +
labs(x = sprintf("Axis1 [%s%%]", round(evals[1], 1)),
y = sprintf("Axis2 [%s%%]", round(evals[2], 1))) +
coord_fixed(sqrt(evals[2] / evals[1])) +
facet_wrap( ~ Season)
# plot_grid(p_ord,p_ord)
print(p_ord4)
# dev.off()
# ggsave(
# paste0(p_ord.file, ".png"),
# p_ord,
# device = "png",
# width = 10,
# height = 6
# )
# gz(paste0(p_ord.file, ".svg"), paste0(p_ord.file, ".svgz"))Test differences between samples on the phylum level
Taxa_tests_phylum1 <- STAMPR(Ps_obj_filt_median, vars2test = "Sample.type", sig_pairs = sig_pairs1, outputfile = paste0(Proj_name, "_Sample.type"))
plotSTAMPR(Taxa_tests_phylum1, pair = "Green tea vs Rooibos")
plotSTAMPR(Taxa_tests_phylum1, pair = "Green tea vs Soil")
plotSTAMPR(Taxa_tests_phylum1, pair = "Rooibos vs Soil")
Taxa_tests_order1 <- STAMPR(Ps_obj_filt_median, vars2test = "Sample.type", rank = "Order", sig_pairs = sig_pairs1, outputfile = paste0(Proj_name, "_Sample.type"))
plotSTAMPR(Taxa_tests_order1, pair = "Green tea vs Rooibos", tax_level = "Order")
plotSTAMPR(Taxa_tests_order1, pair = "Green tea vs Soil", tax_level = "Order")
plotSTAMPR(Taxa_tests_order1, pair = "Rooibos vs Soil", tax_level = "Order")
Taxa_tests_phylum2 <- STAMPR(Ps_obj_filt_median, vars2test = "Season", sig_pairs = sig_pairs2, outputfile = paste0(Proj_name, "_Season"))
plotSTAMPR(Taxa_tests_phylum2, pair = "Winter vs Summer")
Taxa_tests_order2 <- STAMPR(Ps_obj_filt_median, vars2test = "Season", rank = "Order", sig_pairs = sig_pairs2, outputfile = paste0(Proj_name, "_Season"))
plotSTAMPR(Taxa_tests_order2, pair = "Winter vs Summer", tax_level = "Order")
Taxa_tests_phylum4 <- STAMPR(Ps_obj_filt_median_tea, vars2test = "Season", sig_pairs = sig_pairs4, outputfile = paste0(Proj_name, "Teabags_Season"))
plotSTAMPR(Taxa_tests_phylum4, pair = "Winter vs Summer")
Taxa_tests_order4 <- STAMPR(Ps_obj_filt_median_tea, vars2test = "Season", rank = "Order", sig_pairs = sig_pairs4, outputfile = paste0(Proj_name, "Teabags_Season"))
plotSTAMPR(Taxa_tests_order4, pair = "Winter vs Summer", tax_level = "Order")
Differential abundance models
Tag rare phyla (for plotting purposes only)
Ps_obj_filt_media_glom <- tax_glom(Ps_obj_filt_median,
"Phylum",
NArm = TRUE) # glomerate to the phylum level
Ps_obj_filt_media_glom_rel <- transform_sample_counts(Ps_obj_filt_media_glom, function(x) x / sum(x)) # transform to rel. ab.
Ps_obj_filt_media_glom_rel_DF <- psmelt(Ps_obj_filt_media_glom_rel) # generate a df
Ps_obj_filt_media_glom_rel_DF$Phylum %<>% as.character() # factor to char
# group dataframe by Phylum, calculate median rel. abundance
Ps_obj_filt_media_glom_rel_DF %>%
group_by(Phylum) %>%
summarise(median = median(Abundance)) ->
medians
# find Phyla whose median rel. abund. is less than 0.5%
Rare_phyla <- medians[medians$median <= 0.005, ]$Phylum
# change their name to "Rare"
Ps_obj_filt_media_glom_rel_DF[Ps_obj_filt_media_glom_rel_DF$Phylum %in% Rare_phyla, ]$Phylum <- 'Rare'
# re-group
Ps_obj_filt_media_glom_rel_DF %>%
group_by(Phylum) %>%
summarise(Abundance = sum(Abundance)) %>%
arrange(desc(Abundance)) -> Taxa_rankDetect differentially abundant OTUs using ALDEx2 (Fernandes et al. 2013)
Sample type differences
significance = 0.05
# run full model
data2test <- t(otu_table(Ps_obj_filt_median))
# comparison <- as.character(get_variable(Ps_obj_filt_median, "Sample.type"))
#ALDEx_full <- aldex.clr(data2test, comparison, mc.samples = 128, denom = "iqlr", verbose = TRUE, useMC = TRUE) # iqlr for slight assymetry in composition
#ALDEx_full_glm <- aldex.glm(ALDEx_full, comparison, useMC = TRUE) # for more than two conditions
#sig_taxa <- rownames(ALDEx_full_glm)[ALDEx_full_glm$glm.eBH < 0.05] # save names of taxa that are significant under the full model
#write.csv(sig_taxa, file = "Aldex_full_significant_taxa.csv")
# Pairwise comparisons
#
ALDEx_comparisons <- list()
ALDEx_comparisons$Comparisons <- sig_pairs1
# Ps_obj_filt_median_subset <- prune_taxa(names(sort(taxa_sums(Ps_obj_filt_median), TRUE)[1:100]), Ps_obj_filt_median)
for (j in seq(1, length(ALDEx_comparisons$Comparisons))) {
# print(j)
ALDEx_comparisons$Comparisons[j] %>%
str_split(., " vs ", simplify = FALSE) %>%
unlist() ->
comparison_string
Ps_obj_filt_median %>%
subset_samples(Sample.type == comparison_string[1] | Sample.type == comparison_string[2]) ->
# subset_samples(Treatment == comparison_string[2] | Treatment == comparison_string[4]) %>%
# subset_samples(Year == "2016" | Year == "2017" | Year == "2018") ->
Ps_obj_filt_median_pairwise
# Remove species with prevalence < 10%
Ps_obj_filt_median_pairwise_s <- DropRareSpecies(Ps_obj = Ps_obj_filt_median_pairwise, prevalence = 0.1)
sample_data(Ps_obj_filt_median_pairwise_s)$Sample.type %<>% fct_relevel(., str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[1])
# make Joint.sample.name for matching OTUs between compared samples (for GGPlotTopOTUs)
suppressWarnings(
sample_data(Ps_obj_filt_median_pairwise_s) %<>%
as(., "data.frame") %>%
rownames_to_column() %>%
mutate_if(is.factor, as.character) %>%
mutate(Joint.sample.name = paste0(.$Sample.type, "_", .$Field, "_", .$Season, "_", .$Replicate)) %>%
mutate_at("Joint.sample.name", ~str_replace_all(., paste0(unique(comparison_string), collapse = "|"), "")) %>% # remove the levels participating in the comparison from the names
mutate_at("Joint.sample.name", ~str_replace_all(., "^_", "")) %>% # remove first "_"
column_to_rownames()
)
ALDEx2plot_pairwise <- CalcALDEx(
physeq_obj = Ps_obj_filt_median_pairwise_s,
vars2test = "Sample.type",
rare_phyla = Rare_phyla,
sig_level = significance,
LFC = 0
)
ALDEx_comparisons$Results[[j]] <- ALDEx2plot_pairwise # store results
ALDEx_comparisons$Results[[j]]$Var1 <- str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[1]
ALDEx_comparisons$Results[[j]]$Var2 <- str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[2]
ALDEX_summary <- tibble(Label = c(paste0("⬆", sum( ALDEx_comparisons$Results[[j]]$effect > 0 & ALDEx_comparisons$Results[[j]]$Significance == "Pass"), " ⬇", sum( ALDEx_comparisons$Results[[j]]$effect < 0 & ALDEx_comparisons$Results[[j]]$Significance == "Pass"), " (", nrow( ALDEx_comparisons$Results[[j]]), ")")))
ALDEx2plot_pairwise %>%
filter(Significance == "Pass") %>%
dplyr::select(OTU, baseMean, effect, Phylum, Class, Order, Family, Genus) %>%
arrange(desc(abs(effect))) ->
ALDEx2plot_pairwise_results
# print(ALDEx2plot_pairwise_results %>%
# kable(., digits = c(2), caption = "Significantly different taxa:") %>%
# kable_styling(
# bootstrap_options = c("striped", "hover", "condensed", "responsive"),
# full_width = F
# ))
write.csv(ALDEx2plot_pairwise_results, file = paste0("Aldex", "_", paste0(comparison_string, collapse = "_"), ".csv"))
# Plot ALDEX plot
p1 <- GGPlotALDExTax(ALDEx2plot_pairwise, OTU_labels = FALSE, sig_level = significance) +
# ggtitle(ALDEx_comparisons$Comparisons[j]) +
geom_text(
data = ALDEX_summary,
mapping = aes(x = Inf, y = Inf, label = Label),
hjust = 1.1,
vjust = 1.6
)
p1 <- p1 + labs(title = ALDEx_comparisons$Comparisons[j])
print(p1)
# Plot OTU plots
# GGPlotOTU(Ps_obj_filt_subset_pairwise_s, vars2test = "Spill.Treatment", "Seq_8")
p2 <- GGPlotTopOTUs(
physeq_obj = Ps_obj_filt_median_pairwise_s,
vars2test = "Sample.type",
ALDEx_obj = ALDEx2plot_pairwise,
rank_by = "effect",
Ntop = 12
)
print(p2)
}phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 92 samples ] sample_data() Sample Data: [ 92 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 3042 taxa and 92 samples ] sample_data() Sample Data: [ 92 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 3042 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 3042 tips and 3040 internal nodes ] 
phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 74 samples ] sample_data() Sample Data: [ 74 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 3217 taxa and 74 samples ] sample_data() Sample Data: [ 74 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 3217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 3217 tips and 3215 internal nodes ] 
phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 76 samples ] sample_data() Sample Data: [ 76 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 3411 taxa and 76 samples ] sample_data() Sample Data: [ 76 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 3411 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 3411 tips and 3409 internal nodes ] 
filter(ALDEx_comparisons$Results[[1]], Significance == "Pass")$OTU[common_taxa <- filter(ALDEx_comparisons$Results[[1]], Significance == "Pass")$OTU %in% filter(ALDEx_comparisons$Results[[2]], Significance == "Pass")$OTU][1] "Seq_64" "Seq_475" "Seq_631" "Seq_620" "Seq_120" "Seq_112" "Seq_797" [8] "Seq_169" "Seq_37" "Seq_124" "Seq_170" "Seq_1032" "Seq_1382" "Seq_403" [15] "Seq_1094" "Seq_763" "Seq_644" "Seq_855" "Seq_947" "Seq_922" "Seq_507" [22] "Seq_75" "Seq_1003" "Seq_659" "Seq_116" "Seq_698" "Seq_323" "Seq_553" [29] "Seq_313" "Seq_733" "Seq_152" "Seq_968" "Seq_1165" "Seq_838" "Seq_25"
[36] "Seq_2055" "Seq_964" "Seq_434" "Seq_469" "Seq_436" "Seq_65" "Seq_83"
[43] "Seq_1042" "Seq_800" "Seq_488" "Seq_237" "Seq_446" "Seq_264" "Seq_535" [50] "Seq_101" "Seq_665" "Seq_268" "Seq_192" "Seq_940" "Seq_281" "Seq_1305" [57] "Seq_166" "Seq_364" "Seq_2092" "Seq_847" "Seq_868" "Seq_175" "Seq_76"
[64] "Seq_68" "Seq_316" "Seq_122" "Seq_496" "Seq_23" "Seq_195" "Seq_17"
[71] "Seq_15" "Seq_202" "Seq_348" "Seq_270" "Seq_431" "Seq_154" "Seq_61"
[78] "Seq_1232" "Seq_709" "Seq_1048" "Seq_214" "Seq_358" "Seq_1128" "Seq_1001" [85] "Seq_854" "Seq_159" "Seq_42" "Seq_439" "Seq_419" "Seq_33" "Seq_380" [92] "Seq_1058" "Seq_1172" "Seq_367" "Seq_397" "Seq_126" "Seq_132" "Seq_483" [99] "Seq_52" "Seq_306" "Seq_115" "Seq_420" "Seq_31" "Seq_41" "Seq_251" [106] "Seq_12" "Seq_5" "Seq_56" "Seq_2" "Seq_109" "Seq_534" "Seq_131" [113] "Seq_664" "Seq_452" "Seq_208" "Seq_227" "Seq_1404" "Seq_786" "Seq_287" [120] "Seq_235" "Seq_646" "Seq_1057" "Seq_801" "Seq_223" "Seq_635" "Seq_310" [127] "Seq_489" "Seq_318" "Seq_48" "Seq_501" "Seq_438" "Seq_95" "Seq_1738" [134] "Seq_1756" "Seq_2179" "Seq_645" "Seq_1294" "Seq_984" "Seq_788" "Seq_1470" [141] "Seq_228" "Seq_111" "Seq_372" "Seq_253" "Seq_10" "Seq_58" "Seq_1685" [148] "Seq_1186" "Seq_447" "Seq_1099" "Seq_130" "Seq_11" "Seq_578" "Seq_458" [155] "Seq_595" "Seq_565" "Seq_6" "Seq_840" "Seq_2786"
Seasonal differences Comparing the seasonsal differences in teabag samples only
significance = 0.05
data2test <- t(otu_table(Ps_obj_filt_median_tea))
# Pairwise comparisons
ALDEx_comparisons <- list()
ALDEx_comparisons$Comparisons <- sig_pairs2
# Ps_obj_filt_median_tea_subset <- prune_taxa(names(sort(taxa_sums(Ps_obj_filt_median_tea), TRUE)[1:100]), Ps_obj_filt_median_tea)
for (j in seq(1, length(ALDEx_comparisons$Comparisons))) {
# print(j)
ALDEx_comparisons$Comparisons[j] %>%
str_split(., " vs ", simplify = FALSE) %>%
unlist() ->
comparison_string
Ps_obj_filt_median_tea %>%
subset_samples(Season == comparison_string[1] | Season == comparison_string[2]) ->
# subset_samples(Treatment == comparison_string[2] | Treatment == comparison_string[4]) %>%
# subset_samples(Year == "2016" | Year == "2017" | Year == "2018") ->
Ps_obj_filt_median_tea_pairwise
# Remove species with prevalence < 10%
Ps_obj_filt_median_tea_pairwise_s <- DropRareSpecies(Ps_obj = Ps_obj_filt_median_tea_pairwise, prevalence = 0.1)
sample_data(Ps_obj_filt_median_tea_pairwise_s)$Season %<>% fct_relevel(., str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[1])
# make Joint.sample.name for matching OTUs between compared samples (for GGPlotTopOTUs)
suppressWarnings(
sample_data(Ps_obj_filt_median_tea_pairwise_s) %<>%
as(., "data.frame") %>%
rownames_to_column() %>%
mutate_if(is.factor, as.character) %>%
mutate(Joint.sample.name = paste0(.$Sample.type, "_", .$Field, "_", .$Season, "_", .$Replicate)) %>%
mutate_at("Joint.sample.name", ~str_replace_all(., paste0(unique(comparison_string), collapse = "|"), "")) %>% # remove the levels participating in the comparison from the names
mutate_at("Joint.sample.name", ~str_replace_all(., "^_", "")) %>% # remove first "_"
mutate_at("Joint.sample.name", ~str_replace_all(., "__", "_")) %>% # remove double "_"
column_to_rownames()
)
ALDEx2plot_pairwise <- CalcALDEx(
physeq_obj = Ps_obj_filt_median_tea_pairwise_s,
vars2test = "Season",
rare_phyla = Rare_phyla,
sig_level = significance,
LFC = 0
)
ALDEx_comparisons$Results[[j]] <- ALDEx2plot_pairwise # store results
ALDEx_comparisons$Results[[j]]$Var1 <- str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[1]
ALDEx_comparisons$Results[[j]]$Var2 <- str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[2]
ALDEX_summary <- tibble(Label = c(paste0("⬆", sum( ALDEx_comparisons$Results[[j]]$effect > 0 & ALDEx_comparisons$Results[[j]]$Significance == "Pass"), " ⬇", sum( ALDEx_comparisons$Results[[j]]$effect < 0 & ALDEx_comparisons$Results[[j]]$Significance == "Pass"), " (", nrow( ALDEx_comparisons$Results[[j]]), ")")))
ALDEx2plot_pairwise %>%
filter(Significance == "Pass") %>%
dplyr::select(OTU, baseMean, effect, Phylum, Class, Order, Family, Genus) %>%
arrange(desc(abs(effect))) ->
ALDEx2plot_pairwise_results
# print(ALDEx2plot_pairwise_results %>%
# kable(., digits = c(2), caption = "Significantly different taxa:") %>%
# kable_styling(
# bootstrap_options = c("striped", "hover", "condensed", "responsive"),
# full_width = F
# ))
write.csv(ALDEx2plot_pairwise_results, file = paste0("Aldex", "_", paste0(comparison_string, collapse = "_"), ".csv"))
# Plot ALDEX plot
p1 <- GGPlotALDExTax(ALDEx2plot_pairwise, OTU_labels = FALSE, sig_level = significance) +
# ggtitle(ALDEx_comparisons$Comparisons[j]) +
geom_text(
data = ALDEX_summary,
mapping = aes(x = Inf, y = Inf, label = Label),
hjust = 1.1,
vjust = 1.6
)
p1 <- p1 + labs(title = ALDEx_comparisons$Comparisons[j])
print(p1)
# Plot OTU plots
p2 <- GGPlotTopOTUs(
physeq_obj = Ps_obj_filt_median_tea_pairwise_s,
vars2test = "Season",
ALDEx_obj = ALDEx2plot_pairwise,
rank_by = "effect",
Ntop = 12
)
print(p2)
}phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 47 samples ] sample_data() Sample Data: [ 47 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 2700 taxa and 47 samples ] sample_data() Sample Data: [ 47 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 2700 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 2700 tips and 2698 internal nodes ] 
phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 45 samples ] sample_data() Sample Data: [ 45 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 3875 taxa and 45 samples ] sample_data() Sample Data: [ 45 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 3875 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 3875 tips and 3873 internal nodes ] 
phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 48 samples ] sample_data() Sample Data: [ 48 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 3270 taxa and 48 samples ] sample_data() Sample Data: [ 48 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 3270 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 3270 tips and 3268 internal nodes ] 
phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 44 samples ] sample_data() Sample Data: [ 44 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 2511 taxa and 44 samples ] sample_data() Sample Data: [ 44 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 2511 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 2511 tips and 2509 internal nodes ] 
phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 47 samples ] sample_data() Sample Data: [ 47 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 1818 taxa and 47 samples ] sample_data() Sample Data: [ 47 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 1818 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 1818 tips and 1816 internal nodes ] 
phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 45 samples ] sample_data() Sample Data: [ 45 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 2745 taxa and 45 samples ] sample_data() Sample Data: [ 45 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 2745 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 2745 tips and 2743 internal nodes ] 
filter(ALDEx_comparisons$Results[[1]], Significance == "Pass")$OTU[common_taxa <- filter(ALDEx_comparisons$Results[[1]], Significance == "Pass")$OTU %in% filter(ALDEx_comparisons$Results[[2]], Significance == "Pass")$OTU][1] "Seq_511" "Seq_35" "Seq_99" "Seq_119" "Seq_540" "Seq_87" "Seq_477" [8] "Seq_1578" "Seq_169" "Seq_371" "Seq_216" "Seq_983" "Seq_170" "Seq_1351" [15] "Seq_526" "Seq_536" "Seq_644" "Seq_2982" "Seq_1503" "Seq_1904" "Seq_2209" [22] "Seq_222" "Seq_193" "Seq_2377" "Seq_2325" "Seq_552" "Seq_1633" "Seq_321" [29] "Seq_1632" "Seq_550" "Seq_796" "Seq_297" "Seq_1548" "Seq_1581" "Seq_1866" [36] "Seq_1658" "Seq_265" "Seq_965" "Seq_1569" "Seq_151" "Seq_2024" "Seq_2469" [43] "Seq_276" "Seq_1539" "Seq_2269" "Seq_667" "Seq_381" "Seq_340" "Seq_584" [50] "Seq_211" "Seq_365" "Seq_661" "Seq_734" "Seq_615" "Seq_416" "Seq_81"
[57] "Seq_153" "Seq_1669" "Seq_1264" "Seq_1432" "Seq_613" "Seq_1353" "Seq_217" [64] "Seq_1911" "Seq_66" "Seq_791" "Seq_1895" "Seq_780" "Seq_201" "Seq_527" [71] "Seq_267" "Seq_47" "Seq_94" "Seq_9" "Seq_275" "Seq_424" "Seq_741" [78] "Seq_1731" "Seq_606" "Seq_113" "Seq_55" "Seq_16" "Seq_134" "Seq_2503" [85] "Seq_1173" "Seq_29" "Seq_637" "Seq_18" "Seq_27" "Seq_232" "Seq_101" [92] "Seq_608" "Seq_164" "Seq_1326" "Seq_1025" "Seq_1704" "Seq_481" "Seq_149" [99] "Seq_619" "Seq_1118" "Seq_1525" "Seq_486" "Seq_96" "Seq_1024" "Seq_50"
[106] "Seq_286" "Seq_60" "Seq_618" "Seq_364" "Seq_647" "Seq_531" "Seq_790" [113] "Seq_1267" "Seq_2041" "Seq_2228" "Seq_1961" "Seq_1197" "Seq_2482" "Seq_187" [120] "Seq_336" "Seq_979" "Seq_70" "Seq_285" "Seq_175" "Seq_991" "Seq_2045" [127] "Seq_713" "Seq_1108" "Seq_927" "Seq_625" "Seq_1255" "Seq_263" "Seq_582" [134] "Seq_249" "Seq_108" "Seq_832" "Seq_366" "Seq_329" "Seq_215" "Seq_376" [141] "Seq_464" "Seq_567" "Seq_32" "Seq_442" "Seq_1320" "Seq_1471" "Seq_764" [148] "Seq_851" "Seq_1732" "Seq_385" "Seq_257" "Seq_161" "Seq_341" "Seq_107" [155] "Seq_139" "Seq_1071" "Seq_17" "Seq_15" "Seq_54" "Seq_30" "Seq_82"
[162] "Seq_431" "Seq_74" "Seq_1354" "Seq_1579" "Seq_157" "Seq_330" "Seq_71"
[169] "Seq_98" "Seq_459" "Seq_7" "Seq_1613" "Seq_441" "Seq_537" "Seq_1750" [176] "Seq_704" "Seq_752" "Seq_465" "Seq_444" "Seq_312" "Seq_466" "Seq_669" [183] "Seq_224" "Seq_214" "Seq_358" "Seq_600" "Seq_627" "Seq_427" "Seq_180" [190] "Seq_351" "Seq_200" "Seq_210" "Seq_142" "Seq_1903" "Seq_273" "Seq_1226" [197] "Seq_1387" "Seq_127" "Seq_255" "Seq_57" "Seq_191" "Seq_1286" "Seq_199" [204] "Seq_1477" "Seq_13" "Seq_1681" "Seq_42" "Seq_724" "Seq_229" "Seq_33"
[211] "Seq_598" "Seq_414" "Seq_296" "Seq_2011" "Seq_2170" "Seq_2004" "Seq_357" [218] "Seq_189" "Seq_91" "Seq_1113" "Seq_413" "Seq_344" "Seq_575" "Seq_1980" [225] "Seq_360" "Seq_59" "Seq_828" "Seq_2014" "Seq_242" "Seq_2856" "Seq_793" [232] "Seq_31" "Seq_41" "Seq_168" "Seq_230" "Seq_293" "Seq_1053" "Seq_877" [239] "Seq_36" "Seq_39" "Seq_62" "Seq_93" "Seq_205" "Seq_394" "Seq_103" [246] "Seq_73" "Seq_162" "Seq_1997" "Seq_684" "Seq_546" "Seq_530" "Seq_1020" [253] "Seq_1175" "Seq_1257" "Seq_693" "Seq_484" "Seq_642" "Seq_1473" "Seq_188" [260] "Seq_1612" "Seq_288" "Seq_208" "Seq_646" "Seq_2110" "Seq_271" "Seq_184" [267] "Seq_1317" "Seq_318" "Seq_136" "Seq_48" "Seq_146" "Seq_95" "Seq_523" [274] "Seq_158" "Seq_1134" "Seq_811" "Seq_90" "Seq_156" "Seq_308" "Seq_290" [281] "Seq_788" "Seq_2983" "Seq_110" "Seq_1915" "Seq_58" "Seq_1360" "Seq_1050" [288] "Seq_1090" "Seq_1135" "Seq_471" "Seq_369" "Seq_1033" "Seq_11" "Seq_458" [295] "Seq_595" "Seq_2638" "Seq_114"
sessioninfo::session_info() %>%
details::details(
summary = 'Current session info',
open = TRUE
) Current session info
─ Session info ─────────────────────────────────────────────────────────────────────────
setting value
version R version 4.0.3 (2020-10-10)
os Ubuntu 18.04.5 LTS
system x86_64, linux-gnu
ui X11
language en_GB
collate en_DK.UTF-8
ctype en_DK.UTF-8
tz Europe/Prague
date 2021-05-21
─ Packages ─────────────────────────────────────────────────────────────────────────────
package * version date lib source
ade4 * 1.7-16 2020-10-28 [1] CRAN (R 4.0.2)
affy 1.68.0 2020-10-27 [1] Bioconductor
affyio 1.60.0 2020-10-27 [1] Bioconductor
ALDEx2 * 1.22.0 2020-10-27 [1] Bioconductor
ape 5.5 2021-04-25 [1] CRAN (R 4.0.3)
assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.0.2)
backports 1.2.1 2020-12-09 [1] CRAN (R 4.0.2)
Biobase * 2.50.0 2020-10-27 [1] Bioconductor
BiocGenerics * 0.36.1 2021-04-16 [1] Bioconductor
BiocManager 1.30.15 2021-05-11 [1] CRAN (R 4.0.3)
BiocParallel 1.24.1 2020-11-06 [1] Bioconductor
biomformat 1.18.0 2020-10-27 [1] Bioconductor
Biostrings 2.58.0 2020-10-27 [1] Bioconductor
bitops 1.0-7 2021-04-24 [1] CRAN (R 4.0.3)
broom * 0.7.6 2021-04-05 [1] CRAN (R 4.0.3)
bslib 0.2.5.1 2021-05-18 [1] CRAN (R 4.0.3)
cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.0.2)
cli 2.5.0 2021-04-26 [1] CRAN (R 4.0.3)
clipr 0.7.1 2020-10-08 [1] CRAN (R 4.0.2)
cluster 2.1.2 2021-04-17 [1] CRAN (R 4.0.3)
codetools 0.2-18 2020-11-04 [1] CRAN (R 4.0.2)
colorspace 2.0-1 2021-05-04 [1] CRAN (R 4.0.3)
cowplot * 1.1.1 2020-12-30 [1] CRAN (R 4.0.2)
crayon 1.4.1 2021-02-08 [1] CRAN (R 4.0.3)
data.table 1.14.0 2021-02-21 [1] CRAN (R 4.0.3)
DBI 1.1.1 2021-01-15 [1] CRAN (R 4.0.3)
dbplyr 2.1.1 2021-04-06 [1] CRAN (R 4.0.3)
DelayedArray 0.16.3 2021-03-24 [1] Bioconductor
desc 1.3.0 2021-03-05 [1] CRAN (R 4.0.3)
details 0.2.1 2020-01-12 [1] CRAN (R 4.0.2)
digest 0.6.27 2020-10-24 [1] CRAN (R 4.0.2)
doParallel * 1.0.16 2020-10-16 [1] CRAN (R 4.0.2)
dplyr * 1.0.6 2021-05-05 [1] CRAN (R 4.0.3)
ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.0.3)
evaluate 0.14 2019-05-28 [1] CRAN (R 4.0.2)
extrafont * 0.17 2014-12-08 [1] CRAN (R 4.0.2)
extrafontdb 1.0 2012-06-11 [1] CRAN (R 4.0.2)
fansi 0.4.2 2021-01-15 [1] CRAN (R 4.0.3)
farver 2.1.0 2021-02-28 [1] CRAN (R 4.0.3)
forcats * 0.5.1 2021-01-27 [1] CRAN (R 4.0.3)
foreach * 1.5.1 2020-10-15 [1] CRAN (R 4.0.2)
fs 1.5.0 2020-07-31 [1] CRAN (R 4.0.2)
generics 0.1.0 2020-10-31 [1] CRAN (R 4.0.2)
GenomeInfoDb 1.26.7 2021-04-08 [1] Bioconductor
GenomeInfoDbData 1.2.4 2021-05-05 [1] Bioconductor
GenomicRanges 1.42.0 2020-10-27 [1] Bioconductor
ggplot2 * 3.3.3 2020-12-30 [1] CRAN (R 4.0.2)
ggpomological * 0.1.2 2020-08-13 [1] Github (gadenbuie/ggpomological@69f3815)
ggrepel * 0.9.1 2021-01-15 [1] CRAN (R 4.0.3)
ggsci * 2.9 2018-05-14 [1] CRAN (R 4.0.2)
ggthemes * 4.2.4 2021-01-20 [1] CRAN (R 4.0.3)
glue 1.4.2 2020-08-27 [1] CRAN (R 4.0.2)
gridExtra * 2.3 2017-09-09 [1] CRAN (R 4.0.2)
gtable 0.3.0 2019-03-25 [1] CRAN (R 4.0.2)
haven 2.4.1 2021-04-23 [1] CRAN (R 4.0.3)
hexbin 1.28.2 2021-01-08 [1] CRAN (R 4.0.2)
highr 0.9 2021-04-16 [1] CRAN (R 4.0.3)
hms 1.1.0 2021-05-17 [1] CRAN (R 4.0.3)
htmltools 0.5.1.1 2021-01-22 [1] CRAN (R 4.0.3)
httr 1.4.2 2020-07-20 [1] CRAN (R 4.0.2)
igraph 1.2.6 2020-10-06 [1] CRAN (R 4.0.2)
IRanges 2.24.1 2020-12-12 [1] Bioconductor
iterators * 1.0.13 2020-10-15 [1] CRAN (R 4.0.2)
jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.0.3)
jsonlite 1.7.2 2020-12-09 [1] CRAN (R 4.0.2)
kableExtra * 1.3.4 2021-02-20 [1] CRAN (R 4.0.3)
knitr * 1.33 2021-04-24 [1] CRAN (R 4.0.3)
labeling 0.4.2 2020-10-20 [1] CRAN (R 4.0.2)
lattice * 0.20-44 2021-05-02 [1] CRAN (R 4.0.3)
lifecycle 1.0.0 2021-02-15 [1] CRAN (R 4.0.3)
limma 3.46.0 2020-10-27 [1] Bioconductor
lubridate 1.7.10 2021-02-26 [1] CRAN (R 4.0.3)
magrittr * 2.0.1 2020-11-17 [1] CRAN (R 4.0.2)
MASS * 7.3-54 2021-05-03 [1] CRAN (R 4.0.3)
Matrix 1.3-3 2021-05-04 [1] CRAN (R 4.0.3)
MatrixGenerics 1.2.1 2021-01-30 [1] Bioconductor
matrixStats 0.58.0 2021-01-29 [1] CRAN (R 4.0.3)
mgcv 1.8-35 2021-04-18 [1] CRAN (R 4.0.3)
modelr 0.1.8 2020-05-19 [1] CRAN (R 4.0.2)
multtest 2.46.0 2020-10-27 [1] Bioconductor
munsell 0.5.0 2018-06-12 [1] CRAN (R 4.0.2)
NADA * 1.6-1.1 2020-03-22 [1] CRAN (R 4.0.2)
nlme 3.1-152 2021-02-04 [1] CRAN (R 4.0.3)
permute * 0.9-5 2019-03-12 [1] CRAN (R 4.0.2)
phyloseq * 1.34.0 2020-10-27 [1] Bioconductor
pillar 1.6.1 2021-05-16 [1] CRAN (R 4.0.3)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.0.2)
plyr 1.8.6 2020-03-03 [1] CRAN (R 4.0.2)
png 0.1-7 2013-12-03 [1] CRAN (R 4.0.2)
preprocessCore 1.52.1 2021-01-08 [1] Bioconductor
prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.0.2)
progress 1.2.2 2019-05-16 [1] CRAN (R 4.0.2)
ps 1.6.0 2021-02-28 [1] CRAN (R 4.0.3)
purrr * 0.3.4 2020-04-17 [1] CRAN (R 4.0.2)
R6 2.5.0 2020-10-28 [1] CRAN (R 4.0.2)
Rcpp 1.0.6 2021-01-15 [1] CRAN (R 4.0.3)
RCurl 1.98-1.3 2021-03-16 [1] CRAN (R 4.0.3)
readr * 1.4.0 2020-10-05 [1] CRAN (R 4.0.2)
readxl 1.3.1 2019-03-13 [1] CRAN (R 4.0.2)
reprex 2.0.0 2021-04-02 [1] CRAN (R 4.0.3)
reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.0.2)
rhdf5 2.34.0 2020-10-27 [1] Bioconductor
rhdf5filters 1.2.1 2021-05-03 [1] Bioconductor
Rhdf5lib 1.12.1 2021-01-26 [1] Bioconductor
rlang 0.4.11 2021-04-30 [1] CRAN (R 4.0.3)
rmarkdown * 2.8 2021-05-07 [1] CRAN (R 4.0.3)
rprojroot 2.0.2 2020-11-15 [1] CRAN (R 4.0.2)
rstudioapi 0.13 2020-11-12 [1] CRAN (R 4.0.2)
Rttf2pt1 1.3.8 2020-01-10 [1] CRAN (R 4.0.2)
rvest 1.0.0 2021-03-09 [1] CRAN (R 4.0.3)
S4Vectors 0.28.1 2020-12-09 [1] Bioconductor
sass 0.4.0 2021-05-12 [1] CRAN (R 4.0.3)
scales * 1.1.1 2020-05-11 [1] CRAN (R 4.0.2)
sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 4.0.2)
stringi 1.6.2 2021-05-17 [1] CRAN (R 4.0.3)
stringr * 1.4.0 2019-02-10 [1] CRAN (R 4.0.2)
SummarizedExperiment 1.20.0 2020-10-27 [1] Bioconductor
survival * 3.2-11 2021-04-26 [1] CRAN (R 4.0.3)
svglite * 2.0.0 2021-02-20 [1] CRAN (R 4.0.3)
systemfonts 1.0.2 2021-05-11 [1] CRAN (R 4.0.3)
tibble * 3.1.2 2021-05-16 [1] CRAN (R 4.0.3)
tidyr * 1.1.3 2021-03-03 [1] CRAN (R 4.0.3)
tidyselect 1.1.1 2021-04-30 [1] CRAN (R 4.0.3)
tidyverse * 1.3.1 2021-04-15 [1] CRAN (R 4.0.3)
truncnorm * 1.0-8 2018-02-27 [1] CRAN (R 4.0.2)
utf8 1.2.1 2021-03-12 [1] CRAN (R 4.0.3)
vctrs 0.3.8 2021-04-29 [1] CRAN (R 4.0.3)
vegan * 2.5-7 2020-11-28 [1] CRAN (R 4.0.3)
viridisLite 0.4.0 2021-04-13 [1] CRAN (R 4.0.3)
vsn * 3.58.0 2020-10-27 [1] Bioconductor
webshot 0.5.2 2019-11-22 [1] CRAN (R 4.0.2)
withr 2.4.2 2021-04-18 [1] CRAN (R 4.0.3)
xfun 0.23 2021-05-15 [1] CRAN (R 4.0.3)
xml2 1.3.2 2020-04-23 [1] CRAN (R 4.0.2)
XVector 0.30.0 2020-10-27 [1] Bioconductor
yaml 2.2.1 2020-02-01 [1] CRAN (R 4.0.2)
zCompositions * 1.3.4 2020-03-04 [1] CRAN (R 4.0.2)
zlibbioc 1.36.0 2020-10-27 [1] Bioconductor
[1] /home/angel/R/library
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library
References
Fernandes AD, Macklaim JM, Linn TG et al. ANOVA-Like Differential Expression (ALDEx) Analysis for Mixed Population RNA-Seq. PLOS ONE 2013;8:e67019.